program testprog
    use lubricant_class
    use ellipticintegrals_module
    use hertzian_class
    implicit none

    !call general
    !call testlowshearvisc
    !call testlimitingshearstress
    !call testelasticshearmodulus
    !call testhighshearvisc
    !call ellipticintegrals
    call testhertz


contains




subroutine general
    implicit none
    type (lub_t) :: l

    l%rheologybehaviour = viscous
    l%highshearviscositybehaviour = bairviscoplastic
    l%lowshearviscositybehaviour = roelands
    l%lowshearviscosityT0 = 273.15d0 + 40.d0
    l%lowshearviscosity0 = 0.1667
    l%lowshearviscosityalphaZ = 1.0300571753d0
    l%lowshearviscositybetaS = 1.3003671247d0
    l%limitingshearstressdependence = linearpressuredirecttemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 2.d-9
    l%limitingshearstressbeta = 1.d-9

    print *, "l%rheologybehaviour = ", l%rheologybehaviour
    print *, "l%highshearviscositybehaviour = ", l%highshearviscositybehaviour
    print *, "l%lowshearviscositybehaviour = ", l%lowshearviscositybehaviour
    print *, "l%limitingshearstressdependence = ", l%limitingshearstressdependence
    print *, "l%elasticshearmodulusdependence = ", l%elasticshearmodulusdependence
    print *, "l%lowshearviscosityT0 = ", l%lowshearviscosityT0
    print *, "l%lowshearviscosity0 = ", l%lowshearviscosity0
    print *, "l%lowshearviscosityalphaZ = ", l%lowshearviscosityalphaZ
    print *, "l%lowshearviscositybetaS = ", l%lowshearviscositybetaS
    print *, "l%limitingshearstressT0 = ", l%limitingshearstressT0
    print *, "l%limitingshearstress0 = ", l%limitingshearstress0
    print *, "l%limitingshearstressalpha = ", l%limitingshearstressalpha
    print *, "l%limitingshearstressbeta = ", l%limitingshearstressbeta
    print *, "l%elasticshearmodulusT0 = ", l%elasticshearmodulusT0
    print *, "l%elasticshearmodulus0 = ", l%elasticshearmodulus0
    print *, "l%elasticshearmodulusalpha = ", l%elasticshearmodulusalpha
    print *, "l%elasticshearmodulusbeta = ", l%elasticshearmodulusbeta
end subroutine general




subroutine testlowshearvisc
    implicit none
    type (lub_t) :: l

    print *, ""
    print *, "Testing roelands low shear viscosity"
    l%lowshearviscositybehaviour = roelands
    l%lowshearviscosityT0 = 273.15d0 + 40.d0
    l%lowshearviscosity0 = 0.1667d0
    l%lowshearviscosityalphaZ = 1.0300571753d0
    l%lowshearviscositybetaS = 1.3003671247d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 0.1667d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0 + 40.d0, 0.d9), " right answer = ", 0.0262d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.196d9), " right answer = ", 613.65d0
    print *, tangentpressureviscositycoefficient(l, l%lowshearviscosityT0, 0.196d9), " right answer = ", 4.18927895200439369d-8
    print *, tangentpressureviscositycoefficient(l, l%lowshearviscosityT0), " right answer = ", 4.18927895200439369d-8

    print *, ""
    print *, "Testing barusdirecttemperature low shear viscosity"
    l%lowshearviscositybehaviour = barusdirecttemperature
    l%lowshearviscosityT0 = 273.15d0 + 40.d0
    l%lowshearviscosity0 = 0.1667d0
    l%lowshearviscosityalphaZ = 4.18927895200439369d-8
    l%lowshearviscositybetaS = 0.046260909d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 0.1667d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0 + 40.d0, 0.d9), " right answer = ", 0.0262d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.196d9), " right answer = ", 613.65d0
    print *, tangentpressureviscositycoefficient(l, l%lowshearviscosityT0), " right answer = ", 4.18927895200439369d-8

    print *, ""
    print *, "Testing barusinversetemperature low shear viscosity"
    l%lowshearviscositybehaviour = barusinversetemperature
    l%lowshearviscosityT0 = 273.15d0 + 40.d0
    l%lowshearviscosity0 = 0.1667d0
    l%lowshearviscosityalphaZ = 4.18927895200439369d-8
    l%lowshearviscositybetaS = 5115.944d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 0.1667d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0 + 40.d0, 0.d9), " right answer = ", 0.0262d0
    print *, lowshearviscosity(l, l%lowshearviscosityT0, 0.196d9), " right answer = ", 613.65d0
    print *, tangentpressureviscositycoefficient(l, l%lowshearviscosityT0), " right answer = ", 4.18927895200439369d-8
end subroutine testlowshearvisc




subroutine testlimitingshearstress
    implicit none
    type (lub_t) :: l

    print *, ""
    print *, "Testing limiting shear stress in linearpressuredirecttemperature"
    l%limitingshearstressdependence = linearpressuredirecttemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 0.05d0
    l%limitingshearstressbeta = 0.017328679d0
    print *, seconddeviatoricinvariant(1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0), " right answer = ", 3.8729833462074170d0
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, limitingshearstress(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing limiting shear stress in linearpressureinversetemperature"
    l%limitingshearstressdependence = linearpressureinversetemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 0.05d0
    l%limitingshearstressbeta = 1916.359996d0
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, limitingshearstress(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing limiting shear stress in exponentialpressuredirecttemperature"
    l%limitingshearstressdependence = exponentialpressuredirecttemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 3.465735903e-9
    l%limitingshearstressbeta = 0.017328679d0
    print *, seconddeviatoricinvariant(1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0), " right answer = ", 3.8729833462074170d0
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, limitingshearstress(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing limiting shear stress in exponentialpressureinversetemperature"
    l%limitingshearstressdependence = exponentialpressureinversetemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 3.465735903e-9
    l%limitingshearstressbeta = 1916.359996d0
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, limitingshearstress(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, limitingshearstress(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6
end subroutine testlimitingshearstress




subroutine testelasticshearmodulus
    implicit none
    type (lub_t) :: l

    print *, ""
    print *, "Testing elastic shear modulus in linearpressuredirecttemperature"
    l%elasticshearmodulusdependence = linearpressuredirecttemperature
    l%elasticshearmodulusT0 = 273.15d0 + 40.d0
    l%elasticshearmodulus0 = 10.d6
    l%elasticshearmodulusalpha = 0.05d0
    l%elasticshearmodulusbeta = 0.017328679d0
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing elastic shear modulus in linearpressureinversetemperature"
    l%elasticshearmodulusdependence = linearpressureinversetemperature
    l%elasticshearmodulusT0 = 273.15d0 + 40.d0
    l%elasticshearmodulus0 = 10.d6
    l%elasticshearmodulusalpha = 0.05d0
    l%elasticshearmodulusbeta = 1916.359996d0
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing elastic shear modulus in exponentialpressuredirecttemperature"
    l%elasticshearmodulusdependence = exponentialpressuredirecttemperature
    l%elasticshearmodulusT0 = 273.15d0 + 40.d0
    l%elasticshearmodulus0 = 10.d6
    l%elasticshearmodulusalpha = 3.465735903e-9
    l%elasticshearmodulusbeta = 0.017328679d0
    print *, seconddeviatoricinvariant(1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0), " right answer = ", 3.8729833462074170d0
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6

    print *, ""
    print *, "Testing elastic shear modulus in exponentialpressureinversetemperature"
    l%elasticshearmodulusdependence = exponentialpressureinversetemperature
    l%elasticshearmodulusT0 = 273.15d0 + 40.d0
    l%elasticshearmodulus0 = 10.d6
    l%elasticshearmodulusalpha = 3.465735903e-9
    l%elasticshearmodulusbeta = 1916.359996d0
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.d9), " right answer = ", 10.d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0+40.d0, 0.d9), " right answer = ", 5.0d6
    print *, elasticshearmodulus(l, l%lowshearviscosityT0, 0.2d9), " right answer = ", 20.d6
end subroutine testelasticshearmodulus




subroutine testhighshearvisc
    implicit none
    type (lub_t) :: l
    double precision :: temperature, pressure, gammaxz(10) !, gammayz, gammaxy, gammaxx, gammayy, gammazz
    double precision :: lsvisc, taulim, hsvisc(10)
    print *, ""
    print *, "Testing high shear viscosity in highshearviscosity"
    l%rheologybehaviour = viscous
    l%highshearviscositybehaviour = bairviscoplastic
    l%lowshearviscositybehaviour = roelands
    l%lowshearviscosityT0 = 273.15d0 + 40.d0
    l%lowshearviscosity0 = 0.1667d0
    l%lowshearviscosityalphaZ = 1.0300571753d0
    l%lowshearviscositybetaS = 1.3003671247d0
    l%limitingshearstressdependence = linearpressuredirecttemperature
    l%limitingshearstressT0 = 273.15d0 + 40.d0
    l%limitingshearstress0 = 10.d6
    l%limitingshearstressalpha = 0.05d0
    l%limitingshearstressbeta = 0.017328679d0
    temperature = l%lowshearviscosityT0
    pressure =  1.d9
    gammaxz = 10.d0**[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    lsvisc = lowshearviscosity(l, temperature, pressure)
    taulim = limitingshearstress(l, temperature, pressure)
    hsvisc = highshearviscosity(l, temperature, pressure, gammaxz)
    print *, "pressure = ", pressure
    print *, "temperature = ", temperature
    print *, "lowshearviscosity", lsvisc, " right answer = ", 0.1667d0
    print *, "limitingshearstress", taulim, " right answer = ", 10.d6
    print *, "highshearviscosity", hsvisc
    print *, "tau = ", hsvisc * gammaxz
end subroutine testhighshearvisc




subroutine ellipticintegrals
    implicit none
    double precision, parameter :: errtol = 1.d-3
    double precision :: x, y, z, w, k2
    double precision :: a, b, k(10)
    integer :: i, ier
    print *, ""
    print *, "#####################################"
    print *, "           ellipticintegrals"
    print *, "#####################################"
    a = 0.2d0
    b = 0.4d0
    k2 = (a*a) / (b*b)
    x= 0.d0
    y= a*a
    z= b*b
    w = RF(x, y, z, errtol, ier)
    print *, "elliptic integral of the first kind at (",x,",",y,",",z,") = ", w
    print *, "ier = ", ier
    x= 0.d0
    y= a*a/(b*b)
    z= 1.d0
    w = RF(x, y, z, errtol, ier)
    print *, "elliptic integral of the first kind at (",x,",",y,",",z,") = ", w
    print *, "ier = ", ier

    x= 0.d0
    y= a*a
    z= b*b
    w = RD(x, y, z, errtol, ier)
    print *, "elliptic integral of the second kind at (",x,",",y,",",z,") = ", w
    print *, "ier = ", ier
    x= 0.d0
    y= a*a/(b*b)
    z= 1.d0
    w = RD(x, y, z, errtol, ier)
    print *, "elliptic integral of the second kind at (",x,",",y,",",z,") = ", w
    print *, "ier = ", ier

    k = [ 1.000d0, 0.631d0, 0.483d0, 0.400d0, 0.346d0, 0.308d0, 0.279d0, 0.256d0, 0.238d0, 0.222d0 ]
    do i = 1, size(k)
        k2 = k(i)*k(i)
        w = RD(0.d0, k2, 1.d0, errtol, ier) / RD(0.d0, 1.d0, k2, errtol, ier)
        print *, "k = ", k(i), "A/B = ", 1.d0/w
    end do

    print *, "RD(0.d0, 0.5d0, 1.d0)=", RD(0.d0, 0.5d0, 1.d0, errtol, ier)
    print *, "RD(0.d0, 0.5d0-1.d-9, 1.d0)=", RD(0.d0, 0.5d0-1.d-9, 1.d0, errtol, ier)
    print *, "RD(0.d0, 0.5d0+1.d-9, 1.d0)=", RD(0.d0, 0.5d0+1.d-9, 1.d0, errtol, ier)

    do i = 1, 100
        k2 = i * 1.d0 / 100
        print *, "k2 =", k2, "RD(0.d0, 1.d0, k2, ier)/RD(0.d0, k2, 1.d0, ier) =", &
                 RD(0.d0, 1.d0, k2, errtol, ier)/RD(0.d0, k2, 1.d0, errtol, ier)
    end do


    do i = 1, 100
        k2 = 1.d0 * i
        w = contactellipticity(k2)
        print "(F6.0 F6.3)", k2, w
    end do
    print *, "RD(0.d0,1.d0,1.d0,ier) =", RD(0.d0, 1.d0, 1.d0, errtol, ier)
end subroutine ellipticintegrals




subroutine testhertz
    implicit none
    type (hertzian_t) :: c
    print *, ""
    print *, ""
    print *, "CIRCULAR CONTACT"
    print *, "Calculated by elliptic contact code"
    c = newhertzianobj()
    c%contactshape = ellipticalcontact
    c%young1 = 210.d9
    c%poiss1 = 0.3d0
    c%young2 = 210.d9
    c%poiss2 = 0.3d0
    c%Fn = 6.d0
    c%rx1 = 25.4d-3/4d0
    c%ry1 = 25.4d-3/4d0
    c%rx2 = 100.d0
    c%ry2 = 100.d0
    c%alpha = 0.d0
    call calchertz(c)
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
    print *, ""
    print *, "Calculated by spherical contact code"
    c = newhertzianobj()
    c%contactshape = sphericalcontact
    c%young1 = 210.d9
    c%poiss1 = 0.3d0
    c%young2 = 210.d9
    c%poiss2 = 0.3d0
    c%Fn = 6.d0
    c%rx1 = 25.4d-3/4d0
    c%rx2 = 100.d0
    call calchertz(c)
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
    print *, ""
    print *, ""
    print *, "LINE CONTACT"
    print *, "Calculated by cylindrical contact code"
    c = newhertzianobj()
    c%contactshape = cylindricalcontact
    c%young1 = 210.d9
    c%poiss1 = 0.3d0
    c%young2 = 210.d9
    c%poiss2 = 0.3d0
    c%Fn = 6.d0
    c%rx1 = 25.4d-3/4d0
    c%rx2 = 1.d4
    c%width = c%ay
    call calchertz(c)
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "width = ", c%width
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
    print *, ""
    print *, ""
    print *, "ALIGNED ELLIPTICAL CONTACT"
    print *, "Ry smaller than Rx"
    c = newhertzianobj()
    c%contactshape = ellipticalcontact
    c%young1 = 210.d9
    c%poiss1 = 0.3d0
    c%young2 = 105.d9
    c%poiss2 = 0.3d0
    c%Fn = 6.d0
    c%ry1 = 10.d-3
    c%rx1 = 20.d-3
    c%ry2 = 30.d-3
    c%rx2 = 40.d-3
    c%alpha = 0.d0
    call calchertz(c)
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
    print *, ""
    print *, "Rx smaller than Ry (perpendicular case)"
    c = newhertzianobj()
    c%contactshape = ellipticalcontact
    c%young1 = 210.d9
    c%poiss1 = 0.3d0
    c%young2 = 105.d9
    c%poiss2 = 0.3d0
    c%Fn = 6.d0
    c%rx1 = 10.d-3
    c%ry1 = 20.d-3
    c%rx2 = 30.d-3
    c%ry2 = 40.d-3
    c%alpha = 0.d0
    call calchertz(c)
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
    print *, ""
    print *, ""
    print *, "CHECK NEW HERTZIAN OBJECT CREATION"
    print *, "Ry smaller than Rx"
    c = newhertzianobj()
    print *, "Fn = ", c%Fn
    print *, "young1 = ", c%young1
    print *, "young2 = ", c%young2
    print *, "poiss1 = ", c%poiss1
    print *, "poiss2 = ", c%poiss2
    print *, "rx1 = ", c%rx1
    print *, "rx2 = ", c%rx2
    print *, "ry1 = ", c%ry1
    print *, "ry2 = ", c%ry2
    print *, "alpha = ", c%alpha
    print *, "effectivemodulus = ", c%effectivemodulus
    print *, "rx = ", c%rx
    print *, "ry = ", c%ry
    print *, "ax = ", c%ax
    print *, "ay = ", c%ay
    print *, "beta = ", c%beta
    print *, "p0 = ", c%p0
    print *, "delta = ", c%delta
end subroutine testhertz




end program testprog
